function U = poisson_solver_2d(mesh, pde)
A = assemble_matrix_2d(mesh, pde);
b = assemble_vector_2d(mesh, pde);
[A, b] = process_boundary_2d(A, b, mesh, pde);
U = A\b;
end

%% assemble_matrix_2d
function matrix_A = assemble_matrix_2d(mesh, pde)
[gauss_points_ref, gauss_weights_ref] = gauss_lengdre_ref("2d");
matrix_A = sparse(mesh.Nb,mesh.Nb);
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    Jacobi = [V(2,:)-V(1,:); V(3,:)-V(1,:)]';
    J_det = det(Jacobi);
    gauss_weights_loc = abs(J_det) * gauss_weights_ref;
    gauss_points_loc(:,1) = V(1,1) + Jacobi(1,1)*gauss_points_ref(:,1) + Jacobi(1,2)*gauss_points_ref(:,2);
    gauss_points_loc(:,2) = V(1,2) + Jacobi(2,1)*gauss_points_ref(:,1) + Jacobi(2,2)*gauss_points_ref(:,2);
    x_hat = (Jacobi(2,2)*(gauss_points_loc(:,1)-V(1,1)) - Jacobi(1,2)*(gauss_points_loc(:,2)-V(1,2)))/J_det;
    y_hat = (-Jacobi(2,1)*(gauss_points_loc(:,1)-V(1,1)) + Jacobi(1,1)*(gauss_points_loc(:,2)-V(1,2)))/J_det;
    dNdx = basis_function_ref(x_hat, y_hat, mesh.basis_type, "dx");
    dNdy = basis_function_ref(x_hat, y_hat, mesh.basis_type, "dy");
    phix = cellfun(@(x,y) (Jacobi(2,2)*x - Jacobi(2,1)*y)/J_det, dNdx, dNdy, "UniformOutput", false);
    phiy = cellfun(@(x,y) (-Jacobi(1,2)*x + Jacobi(1,1)*y)/J_det, dNdx, dNdy, "UniformOutput", false);
    fc = pde.cor(gauss_points_loc(:,1), gauss_points_loc(:,2));
    for i = 1:mesh.Nlb
        for j = 1:mesh.Nlb
            f = fc .* (phix{j} .* phix{i} + phiy{j} .* phiy{i});
            r = gauss_weights_loc * f;
            matrix_A(mesh.elem(n,j), mesh.elem(n,i)) = matrix_A(mesh.elem(n,j), mesh.elem(n,i)) + r;
        end
    end
end
end

%% assemble_vector_2d
function vector_b = assemble_vector_2d(mesh, pde)
[gauss_points_ref, gauss_weights_ref] = gauss_lengdre_ref("2d");
vector_b = zeros(mesh.Nb,1);
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    Jacobi = [V(2,:)-V(1,:); V(3,:)-V(1,:)]';
    J_det = det(Jacobi);
    gauss_weights_loc = abs(J_det) * gauss_weights_ref;
    gauss_points_loc(:,1) = V(1,1) + Jacobi(1,1)*gauss_points_ref(:,1) + Jacobi(1,2)*gauss_points_ref(:,2);
    gauss_points_loc(:,2) = V(1,2) + Jacobi(2,1)*gauss_points_ref(:,1) + Jacobi(2,2)*gauss_points_ref(:,2);
    x_hat = (Jacobi(2,2)*(gauss_points_loc(:,1)-V(1,1)) - Jacobi(1,2)*(gauss_points_loc(:,2)-V(1,2)))/J_det;
    y_hat = (-Jacobi(2,1)*(gauss_points_loc(:,1)-V(1,1)) + Jacobi(1,1)*(gauss_points_loc(:,2)-V(1,2)))/J_det;
    dN = basis_function_ref(x_hat, y_hat, mesh.basis_type);
    ff = pde.f(gauss_points_loc(:,1), gauss_points_loc(:,2));
    for j = 1:mesh.Nlb
        f =  ff .* dN{j};
        r = gauss_weights_loc * f;
        vector_b(mesh.elem(n,j),1) = vector_b(mesh.elem(n,j),1) + r;
    end
end
end